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ABSTRACT 

A common problem in astronomy is the determination of the time shift between two 
otherwise identical time series of measured flux from a variable source, in short the 
determination of a time lag. Two examples of where this problem occurs are in the 
determination of the Hubble constant from multiple images of gravitationally lensed 
variable quasars and also in the determination of distances to OH/IR stars. It is shown 
here that this problem can be seen as a restricted inversion problem. It therefore can be 
solved using the subtractive optimally localized averages (SOLA) method for inversion 
which has been described elsewhere (Pijpers & Thompson 1992, 1994 ; Pijpers & 
Wanders 1994). 

Key words: methods : data analysis - gravitational lensing ~ stars : distances - 
quasars : individual : QSO 0957+561 



1 INVERSE PROBLEMS 

The problem of determining a time lag between two time se- 
ries of fluxes is very similar to the problem of reverberation 
mapping of active galactic nuclei. The difi'erence lies primar- 
ily in what is known about the so-called transfer function 
(cf. Blandford & McKee 1982 ; Peterson 1993). Essentially 
the problem of reverberation mapping comes from a view 
of an AGN as gas clouds surrounding a variable continuum 
source. The gas clouds re-emit the radiation absorbed from 
the continuum source in spectral lines so that the time lag 
between the variation of the line emission and the contin- 
uum emission is a measure of the difference in path length 
to the observer and hence of the distance from the central 
source of the emitting gas clouds. The transfer function is 
thus related to the distribution of clouds around the nucleus. 
Mathematically, the concept of reverberation mapping leads 
to the integral equation 



Lit) 



dr *(r)C(t 



(1) 



Here L and C are the (velocity integrated) line flux of a 
broad line in the AGN spectrum, and the continuum flux 
respectively. Hence the problem is reduced to the inversion 
of the integral equation to obtain the transfer function If 
in equation (1) for ^(r) a Dirac delta function is substituted, 
*(r) — 5{t — to), equation (1) reduces to L{t) = C{t — to). 
Conversely if two light curves are related by a time delay 
to this is equivalent to equation (1) with a transfer func- 
tion 'l'(r) — S{t — to). In this paper the transfer function is 
therefore assumed to be essentially zero everywhere except 
at the unknown time-lag to : ^'(t) = IS{t — to) where 5 is 
the Dirac delta function. 



Contrary to the problem of reverberation mapping 
where there is an assumed causal relationship between the 
variations in continuum and line flux, in the problem of time 
lag determinations the light curves are not distinguishable 
as driving or responding time series. The inversion method 
should reflect this lack of knowledge in a symmetric treat- 
ment of the two time series. The notation of equation (1) is 
slightly modifled to express this : 



F^''\t) 



F^^Ht) = 



dr /5(r-to)F<'''(t-T) 



(2) 



dr ^5(r + to)F'''(t-r 



The two equations both express that there is a time lag to 
between the two time series F'"' and and therefore 

seem redundant. However both need to be used to ensure a 
symmetric treatment of the two time series in the algorithm. 
J is a constant to allow for a constant non-unity ratio be- 
tween the two time series. Note that I can have any value 
larger than and that the time-lag to can be either posi- 
tive or negative, since it is not known a-priori which light 
curve is leading and which is lagging. Just as in the applica- 
tion of the subtractive optimally localized averages method 
(SOLA) to reverberation mapping (cf. Pijpers & Wanders, 
1994 : hereafter PW) the limits of the integration are set to 
a finite value. The reason for this is that for any measured 
time series its total extent is finite and therefore the range 
over which a time lag can be determined is limited to values 
within this range. 

Equations (1) and (2) are idealized in the sense that 
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finite sampling rates of tlie two liglit curves and finite mea- 
surement errors are not yet explicitly accounted for. This 
is done in the following section. Another problem can be 
that one or both of the light curves arc contaminated by a 
foreground or background source. This can be dealt with, 
however, and the procedure is described in the appendix to 
this paper. 

2 A BRIEF DESCRIPTION OF SOLA 

The strategy of the SOLA method in general is to find a set 

of linear coefficients which, when combined with the data, 
produce the value of the unknown convolved function under 
the integral sign for given value(s) of the time lag. In order to 
do this the time series under the integral sign is interpolated. 
To each measurement of the series outside of the integral sign 
is assigned a partial time scries which is a section of the time 
series under the integral sign. As discussed by Pijpers and 
Wanders (PW) this means that to each measurement L{ti) 
the partial time series consisting of the i"* measurement 
C{ti) and all n previous ones, [C(ti-„), C(ti)], is assigned. 
These partial time series [C{ti-n),C(ti)] for all i form the 
set of base functions for the SOLA method when it is applied 
to equation (1). 

These base functions are used to build a localized aver- 
aging kernel using a set of linear coefficients {d} determined 
for this purpose. In other words if Pj^"'' denotes the (interpo- 
lated) partial time series [C{ti),C{ti+n)], which is normal- 
ized to have an integral of 1, and Fi = F{ti) corrected for 
the normalization of the P/"', then 

^Cii^("'=/C(r,ro) 

' (3) 
^CiFi = V(ro) 

i 

where K, is an integration kernel that is localized around tq 
and ^{to) is the associated localized average of * around 
To. The details of this procedure can be found in the papers 
describing the SOLA method (Pijpers & Thompson 1992, 
1994 hereafter PTl and PT2) and its application to rever- 
beration mapping (PW). In many respects the same tech- 
nique is followed here. Figure 1 is included to assist in the 
visualization of this procedure, where the entire time series 
is plotted as a function of time delay. The entire time se- 
ries under the integral sign is re-plotted 10 times with an 
arbitrary vertical offset for each of the measurements of the 
time series outside the integral sign. The integration lim- 
its — Tmax and Tmax are shown as the vertical dash-dotted 
lines. It is clear that if Tmax is sot to a value larger than or 
equal to half the length of the time series no base function 
will be defined over the entire integration interval, which 
renders the summation in equation (3) meaningless. Limit- 
ing the range of integration to smaller values, and therefore 
the possible range over which a time lag can be determined, 
has the result that sections of the time scries can be used 
as base functions. These are the sections between the verti- 
cal dash dotted lines in figure 1. The first few and last few 
measurements then have partial time series associated with 
them that still do not cover the entire range [—Tmax, Tmax] 
which means that these must be excluded from the analy- 
sis. The remaining points have associated partial time series 
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Figure 1. Example of two time series witli 11 measurements. The 
entire time scries under the integral sign is re-plotted 10 times 
with an arbitrary vertical offset for each of the measurements 
of the time series outside the integral sign. The horizontal scale 
is time delay which means that later measurements fall further 
to the left. The solid part is the part of the time series that is 
actually measured, the dashed parts fall before the first or after 
the last of the measurements. 

that are defined within the dash-dotted window in figure 1. 
The total range of integration [—Tmax, Tmax] must be strictly 
smaller than the total length of the measured time series. 
Tmax is a free parameter of the method and since to is un- 
known it may well be necessary to explore a range of values 
for Tmax and re-do the inversion for each value. Note that 
one can increase Tmax at the cost of having less base func- 
tions, i.e. decrease the height and increase the width of the 
window, and vice versa. To construct a good approximation 
to the target form for the integration kernel from the base 
functions it is desirable to have as many base functions as 
possible, which implies a small Tmax. However, to obtain a 
reliable estimate of the time lag it is necessary to ensure 
that Tmax is larger that the expected to. 

The time series is interpolated using Savitzky-Golay in- 
terpolation (cf. Press et al., 1992). In a Savitzky-Golay filter 
a polynomial of fixed degree is fitted to a moving window 
with a fixed rmmbcr of measured data points. In this imple- 
mentation of the SOLA algorithm the data points are chosen 
symmetrically around the subinterval for which the interpo- 
lation is done. The base functions in this implementation 
of the SOLA method are subsections of the entire time se- 
ries and are known to within the errors propagated from the 
measured points by the polynomial fitting algorithm. 

Contrary to the application of SOLA to reverberation 
mapping the kernel /C will not be designed to be localized. To 
obtain an estimate of the parameter to which is the position 
of the peak of the very narrow transfer function appropri- 
ate for simple time lags it is much better to determine the 
first moment of the transfer function. It is clear that if the 
transfer function is narrow compared to the sampling rate 
or intrinsic time scale of variability of F one cannot hope 
to reconstruct its shape in an inversion. The width of the 
reconstructed transfer function will in this case be almost 
identical to the resolution given by the sampling of the time 
series (cf PW and Pijpers 1994) rather than the real width. 
For simple time lag determinations this is not important, 
however, since only the single unknown time lag to needs 
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to be determined. To do this most conveniently one should 
combine the series of partial light curves of the time series 
under the integral sign into a kernel that is not the usual 
Gaussian but a linear function of r. This moans minimizing 



for two sets of coefficients {c^^"'} and {c^^'''} : 

— Tmax 

, (la) (la) 171 

+ M0 2^q 'c) 'Eij 

— '''max 

, (16) (16) I-, 

+ Mo2^q 'c) 'Eij 



In both cases the target function T for the averaging kernel 
is linear : 



(4) 



T = T 



(5) 



The factor /io is a free parameter which can be used to adjust 
the relative weighting of the errors in the variance-covariance 
matrix Eij and the approximation of the target form. The 
use of this parameter has been described in the papers PTl 
PT2 and in PW. Its purpose is to balance a matching of the 
target kernel function against magnification of the errors 
which are generally opposing aims. An extra constraint is 
imposed on the coefficients which is that : 



^cf") = 
^c(") = 



(6) 



Imposing this constraint ensures that the integral of the av- 
eraging kernel is equal to just as the integral of the target 
kernel. In this way any influence of the even-order moments 
of the transfer function is eliminated. The result of the min- 
imization of (4) is that 



(7) 



The superscript for the coefficients d identifies which of the 
time series a or 6 is under the integral sign, and therefore 
used to build the kernel, and the 1 signifies that the linear 
target kernel (5) was used. When (2) and (7) are combined 



(8) 



the result is : 

= J dTl6{T - to) J2 cf {ti - r) 

W j dTl6{T-to)T 
= Ito 

W j dTjS{T + to)T 

whore to is the estimated time lag. It is clear that the right- 
hand sides do indeed give estimates to of the position to of 
the delta function. The factor I can be determined inde- 
pendently by using the SOLA method with a kernel that is 
constant over the entire range of integration as described in 
PW ; T = l/2rmax instead of T = r. (Note that the con- 
stant is chosen to obtain normalization with integral 1 on the 
interval [— Tmax, Tmax])- For this target form the associated 
coefficient sets are {c^""'} and {c'"''}. For the coefficients 
^(Oa) ^ ^(06) jj.^ ^j^^ constraint (6) the sum should be equal to 
1 instead of 0. Note that the normalization of the base func- 
tions is assumed to have been carried out before the kernel 
is constructed in the minimization of (4). 

The effect of data-errors in the measurements on the 
left-hand sides of the equality in equations (2) is taken into 
account in the usual way (cf. PW). Since the result is a linear 
combination of the data the resulting error estimate can be 
computed trivially. The primary reason to prefer the method 
described here to determining the shape of the transfer func- 
tion is that the magnification of data errors is expected to 
be much smaller with the method described here. Low or- 
der moments of ^ such as the zero order moment I and 
the first order moment Ito can be determined with a higher 
accuracy because there is a roughly inversely proportional 
relation between the resolution width of the localized kernel 
IC and the magnification of data errors. 

The effect of data errors in the measurements under the 
integral sign is not quite as straightforward as for the errors 
outside of the integral sign. The most important point to 
realize is that the kernel that is constructed will in general 
not match perfectly the target function. It is for this rea- 
son that the third equality in the two equations (8) is only 
approximate. This means that the estimates to must be cor- 
rected for such deviations from the target form, which can 
be done with the help of the constructed averaging kernel. 
In practice this moans finding at which to a delta function 
must be placed to obtain the estimate to given the aver- 
aging kernel that was constructed. This is straightforward 
since it requires only a simple function evaluation with the 
constructed averaging kernel as the function. For a perfect 
match between those two functions to = to- For real data, 
with only a limited number of base functions to work with, 
there is always a small correction e due to the deviation 
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of the constructed kernel from the target form. Except in 
those sections where the kernels themselves are discussed, 
from this point the correction e is implicit in the estimates 
to : 



to 



(a) 



(la) ^(6) 



+ e 



(a) 



(9) 



E 



+ e 



ib) 



For a determination of the time lag that is explicitly sym- 
metric in the treatment of the two time scries, the mean of 
the time lags To from the two alternatives is taken. Half the 
difference between the two alternatives Z should be equal 
to zero to within the same errors as apply to Tq. 



2 



to 



(a) 



to 



(b) 



(10) 



where (9) is used to calculate the to- It is immediately clear 
that a time lag To thus determined is invariant under inter- 
change of the two time series. Z can be used as an additional 
tool to gauge the accuracy with which the lag is determined 
in the algorithm since it should be equal to to within the 
errors. If it is not this can bo an indication of contamination 
of the time series by an extraneous source. Its value, together 
with the difference between the two determinations (from in- 
terchanging the time series) of the magnification I, can be 
used to correct for such contamination as is demonstrated 
in the appendix. 
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Figure 2. The light curves that were used for the simulations, a. 
without noise, b. with 2 % noise. 



3 ARTIFICIAL DATA 

To assess the influence of data errors on the algorithm 
and possible systematic effects the method was tried out on 
a set of artificial data. The light curves with and without 
errors are shown in figure 2. This light curve is a smoothed 
version of a continuum light curve reported by Peterson et 
al. (1994) for the active galaxy NGC 5548. The sohd fine is 
the original light curve, the dashed line is the light curve af- 
ter convolution with a Gaussian with a width of 0.1 d and a 
central position of 11.3 d. These numbers were chosen arbi- 
trarily but the requirement for this method that the transfer 
function be sharply peaked is satisfied. Also the central po- 
sition is chosen not to be commensurate with any sampling 
interval. The irregular sampling intervals range from a min- 
imum of 1 d to several days. The second panel of figure 
2 shows both of the light curves with random noise added 
drawn from a normal distribution with a sigma of 2 % of 
the flux. Apart from using error free data (0 %), and 2 % 
random noise a final case with 5 % random noise added to 
the 'fiuxes' is also considered. 

It should be noted that although it is usual to measure 
noise compared to the measured flux, it is actually more ap- 
propriate to compare the noise with the amplitude of the 
actual variations in the light curves. An arbitrary constant 
multiplication factor between the fluxes of the measured 
time series contains no information about a time lag. In 
the artificial case at hand this means that the actual S/N is 
~ 25 for the 2 % case and ~ 10 for the 5 % case. 



The first reconstruction is done taking Tmax = 55 d 
which is about 1/5 of the total length of the time series, 
the second is done taking Tmax = 34 d which is about 1/8 
of the total length of the time series, finally a Tmax = 28 d 
which is about 1/10 of the total length is used. These three 
choices of Tmax were applied to all three pairs of time series 
contaminated by noise at %, 2%, and 5% of the fluxes Fi 
respectively. In each case Tmax is the length of the short- 
est partial time series \ti,ti^k\ that is longer than the value 
of Tmax chosen as input for the algorithm. Figure 3 shows 
the constructed averaging kernels for Tmax = 28 d. Shown 
in table 1 are To and Z, the final column shows the error 
estimates from the propagated data errors. The error esti- 
mate is 0.0 for the first 3 entries because for these no noise 
is added to the time series. For the cases with added noise Z 
is within 2cr of as would be expected. The only exception 
is the case for Tmax = 55 where there is a 2.5ct departure 
from 0. It is clear from the large values of Z for the choice 
of Tmax = 55 d for the error free data and for the data with 
random noise added that this value of Tmax is still too large 
to obtain a reliable estimate of the time lag, although for the 
cases where noise is added the deduced value of To is within 
Icr of the true value. There is only marginal reduction of the 
error estimates when reducing Tmax from 34 d to 28 d and 
at least in this realization of the 5 % noise the estimate of 
the time lag is less accurate for the Tmax = 28 d case. 

In figure 3 plots of kernels constructed from the artifi- 
cial data of figure 2 are shown. In this figure the estimates 
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Figure 3. The constructed averaging kernels for Tmax = 28 d. Left-hand panels are for the case without noise, middle panels have 2 % 
noise, right-hand panels have 5 % noise. The grey area denotes the extent of the error bar on the constructed kernel for the cases with 
noise. The dash dotted lines denote a 2(7 interval centred at to • Note that in the lower panels where the role of the time series is reversed 
with respect to the upper panels, the axes axe also reversed to facilitate visual comparison of the time lag to -I- e determined for each 
case. 



Table 1. The simulations. Error 
weighting ^o = 0.01 
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''"max 


To 
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c^(To) 
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10.2 
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0.0 


0% 


34.0 


10.99 


0.26 


0.0 


0% 


28.0 


11.27 


0.13 


0.0 


2% 


55.0 
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4.5 
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34.0 


11.49 


0.78 
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2% 


28.0 


11.37 


0.85 


0.75 


5% 


55.0 
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4.4 


5% 


34.0 


11.0 


2.6 


2.3 


5% 


28.0 


10.1 


2.7 


2.0 



to are on the ordinate scale and the values corrected for ker- 
nel deviations (cf. equation (9) ) are on the abscissa. It can 
be seen that the kernel suffers somewhat from 'edge effects' 
whore the target kernel is only poorly matched. The con- 
structed kernel is flatter than proportional near the edges of 
the integration range. The effect this has is that at the edges 
of the interval [-Tmax, Tmax] the uncertainty in the time lag 
determination is larger than indicated by the formal errors. 
Also, if the real time-lag to is smaller than ~ —20 d it will 
be overestimated by to, and if the time-lag to is larger than 
~ 20 d it will be underestimated by to. This result indicates 



that not only should Tmax be chosen significantly less than 
half the total length of the measured time series, but also 
significantly larger than the expected time lag. Prom these 
simulations it is clear that given a total extent of the artifi- 
cial time scries Ttot = 267 d, Tmax should be chosen < 55 d. 
A 'safe' criterion seems to be : 

Tmax < (11) 
O 

Furthermore because of inevitable 'edge effects' in the kernel 
the true time lag to should be no more than 0.7— 0.8x Tmax. 
This implies that if one wishes to measure reliably a time 
lag using this method the total length of the measured time 
series should be at least 10 — 12 times the time lag to that 
is expected. 



4 THE SAMPLING STRATEGY OF THE 
TIME SERIES OF QSO 0957-|-561 

A crucial aspect of this method and most other methods 
in use to determine time delays is the interpolation of the 
time series, which is almost always necessary to estimate 
time delays. Interpolation schemes arc usually very sensitive 
to the sampling of the time series. Before continuing with 
SOLA time delay determinations it is useful therefore to 
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Figure 4. The optical light curves of images A and B of the grav- 
itationally lensed QSO 0957+561. The fluxes for the two images 
are in the same units chosen such that the values are of order 
unity. (Schild & Thomson, 1995) 

examine in more detail the influence of the samphng on the 
SOLA method. 

4.1 The distribution of measurement points in 
time 

The largest homogeneous data set for the two optical images 
of the gravitationally lensed QSO 0957+561 is that reported 
by Schild & Thomson (1995), who also have made available 
a master set which combines data from other sources. The 
time scries for the former are shown in the two panels of 
figure 4. The total extent of these time series is 5347 d and 
so a time lag to of up to ~ 550 d can be determined with 
reasonable accuracy with this method and these time series. 
Rather than use magnitudes the time series are converted to 
(arbitrary) flux units. The reason for this is that if there is 
any contamination present in the time series, due to extra- 
neous factors such as foreground or background sources or 
micro-lensing, this is additive in flux and not in magnitude. 
Such contamination can be corrected for to some extent as 
long as it is additive, as demonstrated in the appendix. The 
flux is calculated from the magnitudes using : 

Fi = 4 10'' oxp (-mi/2.5) (12) 

Schild and Thomson (1995) note that since the more 
recent part of their measured data series has sections with 



sampling rates that are as high as once per day, there should 
be no lack of high frequency signal in the time series with 
which to determine a time lag with high accuracy. This is 
true in principle but it docs depend to some extent on the 
manner in which these sections are distributed within the 
overall time series. To illustrate this point all those measure- 
ment dates were isolated from the time series for which there 
are also measurements on the previous two days and the fol- 
lowing two days. Thus all these points arc the middle point 
of a section of 5 measurements done on consecutive days. In 
the entire time series there are 366 such quintuplets, many 
of which partially overlap because for some periods there are 
even more than 5 consecutive days on which measurements 
are available. In order to enable using the high frequency 
signal in these quintuplets when determining the time lag 
they must overlap with another quintuplet after the time 
series is shifted by that time lag to- If the time separation 
between a given pair of quintuplets is the exact time lag be- 
tween the two quasar images the measured flux in the two 
images should behave exactly the same (in the absence of 
measurement errors), so the quintuplets should be exactly 
the same. If two quintuplets do not overlap after shifting, 
whatever high frequency signal they contain does not have 
a measured counterpart in the other time series and is virtu- 
ally useless. The time separation between the middle points 
of each of these pairs is the time lag for which that pair can 
be used optimally to determine whether it is the actual time 
lag between the two quasar images. The more quintuplets 
overlap for all possible time lags the better the true time 
lag can be determined with the measured series. Of course 
it is quite arbitrary to choose 5 consecutive days. One could 
as well take any other number or take all weeks for which 
there are a certain minimum number of measurements done. 
This will change some of the histograms to be shown in this 
section but not the essentials of the arguments in favour of 
carefully designing sampling strategies. 

With these = 366 quintuplets in the measured time 
series ^N{N — 1) = 66795 distinct pairs of distinct quintu- 
plets can be formed. For each distinct quintuplet pair the 
time difference between the middle points is calculated and 
then binned into 5 day bins thus accounting for the fact 
that even partially overlapping high sampling rate sections 
are useful. The result is depicted £is a histogram of num- 
ber of pairs per time lag bin. It is clear that if there are 
many overlapping quintuplets for a given time separation 
then a time lag in that range can be well determined by 
the high sampling rate quintuplets, and only poorly if there 
are very few overlapping quintuplets. For an ideal sampling 
strategy, without any prior knowledge of the actual time lag, 
the quintuplet pair separations should cover uniformly the 
entire range of time lags of interest. 

In figure 5 the resulting histogram is shown in the bot- 
tom panel. The largest separation between two quintuplets 
in the measured time series is 3514 d but the plotted his- 
togram is restricted to time lags between and 750 d which 
covers the range of interest since the method presented here 
cannot reliably recover time lags outside of this interval for 
the time series measured for QSO 0957+561. For comparison 
in the top and middle panels two artificially designed sam- 
pling strategies are shown. In the top panel the same number 
of 366 quintuplets have been distributed uniformly over the 
time period covering 3514 d. It can be seen from figure 5 
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Figure 5. Lower panel : A histogram of the time separations 
of quintuplets of consecutive measured points in the time series 
of the gravitational lens QSO 0957+561. The upper panel shows 
what this histogram would be if the quintuplets were uniformly 
spaced over the same time span as the largest quintuplet separa- 
tion in the actual time series. The middle panel shows what this 
histogram would be if the quintuplet separation was designed to 
cover this time span in a non-redundant fashion. 

that this results in a bimodal pattern of separations. For 

some time lags no quintuplet overlaps with another and for 
some lags quite a large number overlaps. The few interme- 
diate bins occur because the separation between successive 
pairs is not an integer multiple of the bin width of 5 d. The 
second strategy shown in the middle panel of figure 5 does 
rather better. Here the A'^ = 366 quintuplets are distributed 
over the same time period covering 3514 d by putting them 
at the times Ts(i) obtained by taking : 

PiPN+l 



X 3514 d V 



T..,iV 



(13) 



where the Pi are prime numbers. Of course the resulting 
measurement times do not fall on integer dates, so each mea- 
surement time is rounded to the nearest integer before the 
separations are calculated and binned. This strategy results 
in a much more nearly uniform coverage of the range of pos- 
sible time lags. The actual distribution of quintuplet pair 
separations shown in the bottom panel of figure 5 is clearly 
not uniform. There are 'gaps' around 150 d and between 
200 d and 250 d and between 470 d and 530 d. Annual oc- 
cultation of the source by the Sun causes large gaps in the 
time series which come back every year. This can diminish 
the number of available quintuplet separations between time 



lags of around the time corresponding to the length of that 
gap and around 1 yr minus that length. Similar dips in the 
distribution of quintuplets separations would occur at the 
delays of these same times with integer numbers of years 
added. However one can compensate partially for this effect 
on the distribution of separations by placing the quintuplets 
more often at the separations that currently seem to be ne- 
glected. For a sun constrained measurement gap that lasts 
less than half a year this is always possible, weather and 
telescope scheduling permitting. 

If the true time lag is between 350 d and 450 d then 
the actual samiiling of the time series appears fortunate, 
since in this range there are very many pairs that overlap. 
If the true time lag is between 450 d and 550 d however, 
the number of overlapping quintuplet pairs decreases to well 
below what it would be in the nearly uniform case shown in 
the middle panel of figure 5. Since there are claims in the 
literature (cf. Press et al., 1992a,b) that the true time lag lies 
in this range the sampling of this time series seems rather 
unfortunate, since the time series is not best suited to test 
the longer lags. However, one must keep in mind that this is 
only true in so fax as the very high sampling rate episodes 
are concerned. The overall time series can of course be used 
for time lag determinations. The distribution of the high 
sampling rate sections does demonstrate that the variations 
seen in those sections can in general only with difficulty be 
treated as signal since their counterparts in the time series 
for the other image are in general not very likely to have 
been measured as well. 

One should not conclude from this discussion that 
SOLA or any other method is unable to recover a time de- 
lay if the sampling times are not distributed uniformly or 
according to (13). The simulated light curve of section 3 
did not follow this sampling and it is clear that a time de- 
lay can be recovered. However the most accurate recovery 
for a given number of samplings and a given photometric 
accuracy is obtained when the separations in time between 
each pair of measurement points are designed to be non- 
redundant. Weather conditions will always prevent the use 
of a pre-designed strategy but by monitoring the measure- 
ment point-pair separations while monitoring and schedul- 
ing subsequent measurements to compensate for any dips in 
the distribution function of these separations it should be 
possible to obtain a much more nearly uniform distribution. 
In the following section the differences between the various 
strategies is demonstrated using simulated time series. 

4.2 Comparing the results for difFerent sampling 

strategies 

To demonstrate the effect of sampling strategies three sam- 
pling strategies are used for the same simulated time series. 
For this simulation the time series for image A is interpo- 
lated linearly and resampled uniformly, and also resampled 
according to (13). The second time series was obtained by 
shifting the original by 511 days. Finally on all these six time 
series random numbers drawn from a Gaussian distribution 
with zero mean and a standard deviation of 1% of the flux 
were added to mimic measurement errors. This simulation 
is somewhat more realistic than the first example discussed 
in section 3 because there is more high frequency signal in 
the image A (and B) time series. 
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Figure 6. The constructed kernels for the three different samphng strategies of a simulated time series. The true delay is 299 d. The 
error weighting // = 10~^ in all cases. The Tmax is 711 for the original sampling of the data, 715 for the unform sampling and 730 for 
the prime number sampling. Only the relevant subsection of the kernel is shown in order to demonstrate the differences between kernels. 



An inversion is carried out using the procedures out- 
lined before. In all cases the Savitzky-Golay fitting of the 
time series is done with a window that is three points wide 
and a constant (0-order polynomial) is fitted to these points. 
The error weighting fj, = 10~^ in all cases. The value of Tmax 
is 711 d for the original sampling of the data, 715 d for the 
uniform sampling and 730 d for the sampling according to 
(13). The results are summarized in table 2. The two ar- 
tificial time lags were chosen to demonstrate the difference 
between the accuracy for a time lag that should be well sam- 
pled even with the original sampling of the A image, and a 
time lag that should be more difficult to recover. The con- 
structed linear kernels for the case where the true To = 299 d 
are shown in figure 6. The results in table 2 show that there 
is negligible difference between the uniform and prime num- 
ber sampling results. From the figure 6 it is clear however 
that the prime number strategy produces a more uniform 
kernel. In the kernel for the regularly sampled time series 
the bimodal pattern of quintuplet separations is refiected in 
the regular block structure of the error bars on the kernel. 
An encouraging result is that even for the original sampling 
the true time lag is quite well recovered. In fact for this re- 
alization of the artificial errors the recovery is marginally 
better for the longer time lag. More worrisome is that the 
asymmetry in the results when interchanging the role of the 
time series is quite large for the true sampling. Since the 
asymmetry is used to determine possible contamination of 



the time series by extraneous efi^ects one might mistakenly 
conclude that the time series is contaminated. Subsequently 
correcting for this non-existent contamination will produce 
a time lag that is rather less accurate than the formal errors 
indicate. It is in particular in the differences 61 and Z that 
one can see that a time lag inversion in the region of 500 d 
is more problematic than one in the region of 300 d. 

5 THE TIME LAG FOR QSO 0957-|-516A,B 

5.1 The Schild & Thomson data 

SOLA time lag determinations are carried out using various 
parameter settings and with the time series obtained for 

the gravitational lens QSO 0957-1-561 recently reported by 
Schild and Thomson (1995). Parameters for the inversion 
are : 

'7'max . 

- The number of points in the moving window to inter- 
polate the time series under the integral sign N^-in- 

- The degree of the polynomial to fit to these points A^'poi. 

- The error weighting parameter /j,. 

The parameter settings for 3 cases are summarized in table 
3. The output of the SOLA method would then be the zero 
order moment /, and the asymmetry SI in this when inter- 
changing the role of the time series, and further To and Z. 
Using the method described in the appendix to correct for 
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Table 2. The results from the time lag inver- 
sions for the artificial time series using different 
sampling strategies. The entries in brackets are 
the propagated data errors which always apply 
to the last decimal place of the entry to the left 
of it. The results for two different artificial lags 
are shown. The columns show in order the mean 
of the two determinations of /, and half the dif- 
ference SI between them, and To, and Z. 



case 


I 


SI 


To 


Z 


true 


1.0000 


0.0000 


299 





original 


0.9927 (4) 


0.0005 (4) 


282 (6) 


23 (6) 


uniform 


0.9990 (3) 


0.0003 (3) 


295 (4) 


8 (6) 


prime no. 


0.9988 (3) 


0.0001 (3) 


294 (5) 


8(6) 


true 


1.0000 


0.0000 


511 





original 


0.9940 (4) 


0.0023 (4) 


514 (6) 


66 (6) 


uniform 


0.9990 (3) 


0.0002 (3) 


515 (4) 


7(6) 


prime no. 


0.9990 (3) 


0.0003 (3) 


511 (5) 


9(6) 



contamination, SI and Z are replaced by a relative oflFset 
2/1 = C(") - C**"'// and a relative drift 2/2 = L^"' - L^''^/I. 
The output is thus four numbers : 

- The relative offset yi. 

- The relative drift 1/2, assumed to be linear in time. 

- The relative magnification I. 

- The time-lag To. 

The results for the three cases considered hero arc summa- 
rized in table 4. Inversions with other parameter settings 
have been carried out. The sensitivity of the errors to the 
error weighting is small and no systematic shifts in the re- 
sults have been found. The value of Tmax is quite strongly 
constrained by the time lag itself and by the total length 
of the measured time series. Other values did not yield sig- 
nificantly different results and usually had larger associated 
error estimates. The cases shown here can thus be consid- 
ered representative of the most reliable determinations done 
with this method and these data. 

It is immediately clear from the asymmetries obtained 
in both the values of / and of the time lags to when revers- 
ing the role of the time series that some contamination of 
the time series is present. This asymmetry is much larger 
than in the artificial time series with the same sampling. 
The cause of the contamination can be found in e.g. inac- 
curate subtraction of a foreground or background source, or 
micro-lensing (cf. Schild and Smith 1991, Falco et al. 1991b). 
Following the procedure outlined in the appendix the asym- 
metries are used to obtain an estimate of this contamination. 
As described in the appendix it is impossible to determine 
from this procedure which time series is contaminated to 
what extent, since only a differential contamination can be 
determined. The contamination in each time series is as- 
sumed to be of the form C + L{t — ti). In table 4 the offset 
2/1 = C^"' - C^''^/I and the drift 2/2 = i'"' - L^''^/I. The 
flux units for 2/1 and 2/2 are the same arbitrary units as used 
in equation (12). Note that now the asymmetries 51 and Z 
are used in effect to determine y\ and 2/2 so that after the 
corrections the SI and Z are equal to to within lu. 

Assuming for convenience that this contamination is en- 
tirely contained in the A image the appropriate contribution 
is subtracted from the time series of image A. Adding a sim- 
ilar contribution to image B instead produces identical re- 
sults in the subsequent inversions. The resulting time series 
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Figure 7. The time scries for image A and B of the gravitational 
lens QSO 0957+561 after removing a linear function of time ex- 
pressing the offset and linear drift due to extraneous sources (see 
the appendix). Only the difference in contamination between A 
and B can be determined so the contamination is arbitrarily sub- 
tracted from only the A image time series. 

for case 2 are plotted in figure 7. After this subtraction the 
inversion yields the relative magnification I and the time 

delay To. The uncertainties quoted in table 4 are la er- 
ror bars due to the propagation of the measurement errors 
of the time series. However one should note that the error 
cross- correlation coefficients for these four numbers are all in 
excess of 0.8. In terms of a measure the iso-x^ contours in 
the 4-dimensional solution space form hyper-ellipsoids with 
their major axes not aligned with the coordinate axes. Tak- 
ing this into account the actual uncertainty in for instance 
the time lag is considerably larger, because an equal can 
be achieved by simultaneously changing the other 3 values. 
Furthermore one should note that the contamination is as- 
sumed to be at most linear in time. As demonstrated by 
Schild and Smith (1991), and by Falco et al. (1991b) there 
may well be higher-order terms present due to the effect of 
micro-lensing. While higher-order contamination will behave 
more as noise rather than cause systematic shifts in the re- 
sult, the quoted measurement error is then no longer a good 
measure of the actual noise in the time series. This effect also 
increases the uncertainty of the results. The spread of points 
within quintuplets is roughly a factor of 2 times the quoted 
error bars. If one assumes that the high frequency micro- 



10 F.P. Pijpers 



Table 3. The various parameter settings in the time lag inversions. For case 5 
all points in quintuplets of daily sampling have been merged as described in the 
text. 



case 


Tmax 


N ■ 

' win 


^pol 






yi 


y2 

[xlO-3 yr-i] 


/ 


To 
[d] 


1 


711 


2 


1 


10" 


3 


0.28 ± 0.02 


-5.8 ±0.3 


1.34 ±0.06 


472 ± 73 


2 


711 


3 





10" 


3 


0.29 ±0.01 


-6.8 ±0.2 


1.35 ± 0.03 


425 ± 17 


3 


711 


5 


1 


10" 


3 


0.28 ± 0.02 


-6.8 ±0.2 


1.30 ±0.05 


434 ± 16 


4 


711 


11 


2 


10" 


3 


0.2 ±0.1 


-12 ±3 


1.2 ±0.3 


377 ± 41 


5 


711 


3 





10" 


3 


0.21 ±0.02 


-6.8 ±0.2 


1.21 ±0.05 


469 ± 16 



lensing events aeeount for this excess spread one should in- 
crease the errors by some amount, keeping in mind that 
such micro-lensing events do not follow the same statistical 
distribution as measurement errors. 

The contamination, the relative amplification factor 
and the time lag shown in Table 4 for cases 1 , 2 and 3 demon- 
strate that the inversions are consistent to within the errors. 
This is what would be expected for an ideal case where the 
errors are Gaussian and the distribution in time of the sam- 
pling points is adequate for a time lag determination. To 
assess the influence of possible high-frequency contaminat- 
ing signal from micro-lensing cases 4 and 5 were carried out. 
In case 4 the interpolating window is set to 11 points to 
which a quadratic polynomial is fitted. This should smooth 
out some of the variation within the window but of course 
the detrimental efi^ect of large gaps in the time series are 
exacerbated by a wider window. An alternative approach is 
to isolate the quintuplets of points measured on consecutive 
days and replace them with their mean on the central day of 
the quintuplet. This is also a smoothing operation but now 
carried out only on sections with a high sampling rate. This 
is done in case 5. Both for case 4 and 5 it is important to re- 
alize that this stronger smoothing is not specific. True high 
frequency variations that arc intrinsic to the Icnsed quasar 
are removed together with the micro-lensing variations. 

As a result of the smoothing one expects that the offset 
2/1 becomes smaller in absolute value. The relative magnifi- 
cation factor is affected as well however. This can be under- 
stood by writing the smoothing in terms of an extra con- 
volution with a smoothing function S{t) that is determined 
by the choice of window and fitting polynomial. Instead of 
having a convolution as in equation (1) : 

L{t) = * C{t) (14) 

the convolution is now : 

L{t) = * S-\t)) * {S{t) * C{t)) (15) 

(or the analogous forms for (2) ). Here? S'^^(t) is the inverse 
of the smoothing operation S. That is to say that because 
the fit S{t) * C{t) is used the transfer function ^ is modi- 
fied to ^'(i) * S~'^{t). When this is used in (2) and the steps 
outlined before arc followed to obtain the relative magnifi- 
cation / it is clear that what is actually obtained is either 
IS~^{to) or I/S~^{to). Increasing the smoothing leads to 
an increasing difference between these two values. This in- 
creased asymmetry leads to an overestimate of yi and an / 
that will tend towards 1. In the limit where all variations 
in the time series are smoothed out one cannot distinguish 



between the cases of time scries that arc offset duo to con- 
tamination and of time series that have a non-unity ratio 
between the flux levels, which is the ratio between the aver- 
age flux in either series. The time lag and the drift y2 should 
be less affected by the smoothing as long as the smoothing 
window is symmetric around the point at which the interpo- 
lation is evaluated. The error in y2 and To should increase 
however because the accuracy with which a lag can be de- 
termined depends in part on the intrinsic high frequency 
signal. The results in Table 4 bear out the expectations in 
that not only yi but also I appear reduced. Because the 
2/1 is probably overestimated the reduction of the contam- 
ination by the smoothing is actually larger than the result 
shows because it is partially negated. Considering the errors 
the various effects remain marginal however. It is not clear 
whether the contamination is entirely due to high frequency 
variations. The drift 2/2 appears to be still present and the 
effects of the smoothing should play only a minor role here. 
The time lag for cases 4 and 5 arc consistent with the best 
determinations but the error is clearly higher for case 4. On 
the basis of these results it seems that the best determina- 
tion of the time lag is case 2. An error weighted mean of all 
5 time lag determinations yields To = 441 ±9 d. Considering 
that these five determinations arc not independent because 
the same data is used it seems preferable to take the best 
case To = 425 ± 17 d as final result. 

5.2 Other data sets 

In the literature, there are reported a number of attempts 
at determining the time lag from optical (Schild & Cholfin 
1986, Schild 1990, Pelt et al. (1994, 1996) at R band and B 
band by Florentin-Nielsen (1984), Vanderriest ct al. (1989), 
Press et al. (1992a) Beskin & Onkyanskij (1992) as well 
as radio data (Roberts et al 1992. Lehar et al 1989, 1992, 
Press et al 1992b) . There is even an attempt using UV data 
(Gondhalekar et al. 1986). The B band monitoring produces 
data that is similar to that used here although generally the 
time series are shorter. Considering only the B band data 
there appear to be emerging two mutually exclusive time 
lags, that are obtained from the various analyses. Press et 
al. analysed a shorter optical time series for QSO 0957+561 
(1992a) and also a radio time series (1992b) using a rigorous 
statistical method developed by them (1992a). They obtain 
an estimate for the time delay of 5361^2 d as 95 % confidence 
interval from the optical data and their determination from 
the radio data is consistent with this value. An analysis by 
Pelt and coworkers (1994) of the same data yields 415 ± 
32 d and their analysis of the radio data is consistent with 
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their value for the optical time series. The determination 
reported here appears at face value to be consistent with the 
determinations of Pelt and coworkers (1994, 1996), and not 
with that of Press ot al. Considering the fact that the SOLA 
method can only correct for offsets and linear drifts and also 
that the correlation coefficient between the yi, j/2, I, and To 
are all large means that it is not justified to exclude the 
estimates of Press et al. out of hand. However Thomson and 
Schild (1996) argue that the presence of variations due to 
micro-lensing invalidates some of the assumptions made by 
Press et al. (1992a, b) concerning the statistical properties of 
the time series. Thomson and Schild (1996) discuss how the 
Press et al. method can be corrected for this effect and then 
obtain a result that is consistent with the shorter time lag 
found previously by them. Since the study presented here 
also shows evidence for contamination of the time series, 
although no specific cause can be identified from the SOLA 
method itself, there is some support for their hypothesis. 

To compare with previous determinations the SOLA 
method is also applied to the same data (Vanderriest et al., 
1989) that the method of Press et al. was applied to. This 
data set spans a period of a bit over 2900 days which is on the 
short side for a reliable estimate of a time lag using SOLA, 
as argued in section 3. Furthermore the errors are larger for 
this data set than for the Schild and Thomson data which 
leads to a larger uncertainty in the time delay as well. The 
same parameter settings are used for the SOLA algorithm as 
in case 3 for the Schild and Thomson data except that Tmax 
is reduced to 601 d. The resulting relative magnification is 
1.0±0.14 and the time delay is 429±49 d. The asymmetries 
51 and Z are on the level of the la errors. 

Finally Schild and Thomson (1995) present a master set 
of data from various sources, including their own. This data 
set spans about the same time period that Schild & Thom- 
son's own data do, but the coverage is somewhat better. The 
SOLA parameter settings are as in case 2 in the previous sec- 
tion. The resulting relative magnification / = 1.30 ± 0.03. 
and the time delay To = 496±12 d. The contamination offset 
is yi = 0.26 ±0.01 and the drift j/2 = -6.4 ±0.2 x 10"^. Al- 
though this set contains more measurements than the others 
considered here one should realize that these data are not 
homogeneous in quality. If the contamination in the time 
series is due in part to imperfect data reduction procedures, 
which differ between the various observers, then it is a very 
poor approximation to assume that the contamination can 
be described in terms of a simple offset and drift. Further- 
more the merging process itself of combining these different 
parts of the time series can introduce extra systematic er- 
rors. Taking this into account it seems unsurprising to find 
a difference of some 5 — 7a between this result and the result 
when using only the data of Schild and Thomson (1995). 

6 OBTAINING Ho FROM THE TIME LAG 

From a careful analysis of the images of QSO 0957+561A, 
B combined with a tracing of the light deflection through a 
model potential to fit the positions and magnifications of the 
images it is possible to determine the gravitational potential 
of the lensing object. With this potential it is possible to 
determine the Hubble constant as long as the Icnsed source 
is variable. This is essentially a consequence of the light of 
the two images having traversed a different path through 



this gravitational potential. By measuring the time delay 
for a signal propagating at the speed of light an apparent 
separation of the images can be related to a physical distance 
which together with a measured redshift yields an estimate 
of the Hubble constant. Further details can be found in the 
monograph 'Gravitational Lenses' (Scheider et al. 1994) 

Prom modeling of the lensing system Falco et al. (1991a) 
quote for the value of the Hubble constant : 

Ho = (90 ± 10) ( ^ X ( ^] ' (16) 

^ '^V390 kms-iy \lyrj ^ ' 

The units are km s^^ Mpc^^, and ay is the velocity- 
dispersion of the lensing galaxy, tab is the time lag between 
the time series of the images. More recent modeling by Gro- 
gin and Narayan (1996) yields the result : 

Using the value av = 303 ± 50 obtained by Rhee (1991), a 
time delay of 425 ± 17 d, and using formula (16) yields : 

i?o = 47 ± 16 km s"^ Mpc"^ (18) 

Using instead (17) yields : 

Ho = 69 ± 24 km s"^ Mpc"^ (19) 

The largest contribution to the error in these determinations 
is due to the velocity dispersion. Grogin and Narayan (1996) 
also quote an equation similar to (17) for an upper limit to 
the Hubble constant which is independent of the velocity 
dispersion : 

Ho = {9>2.^tll){l-K)(^^^^ kms-^Mpc-i (20) 

where k is the convergence of the quadratic potential with 
which the cluster surrounding the lensing galaxy is approxi- 
mated. Since k must be positive, using the time delay yields 
an upper limit : 

Ho < 78 ± 7 km s"^ Mpc"^ (21) 

As Grogin and Narayan (1996) discuss the of their model 
for the gravitational potential is somewhat lower than that 
of Falco et al. (1991a) when using the most recent VLBI 
data as constraints for the model parameters. For this reason 
equations (20) and (17) are probably more reliable than the 
older values of Falco et al. (1991a). 

7 CONCLUSIONS 

The SOLA method for inversion is modified to determine 
the time lag between two time series, assuming that the 
width of the transfer function between the two is small com- 
pared to the typical sampling time interval. The method is 
demonstrated to perform well on artificial data. The method 
is then applied to the measured time series of the gravita- 
tional lens QSO 0957-1-561. Considering the uncertainties in 
the de-trending the best estimate obtained here gives a value 
for the time delay of 425 ± 17 d. This leads to a best value 
for the Hubble constant of Hq = 69 ± 24 km s"^ Mpc"^ 

Considering the primary sources of uncertainty in the 
determination of the Hubble constant from the time delay 
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for QSO 0957+516 it appears highly desirable to obtain a 
more accurate determination for the velocity dispersion of 
the lensing system. The determination of the time lag it- 
self could be improved upon substantially if an independent 
means can be found to determine quantitatively the effects 
that contaminate the time series such as micro-lensing but 
also the effects of merging data from various sources and 
observers. The SOLA time lag determinations from Schild 
& Thomson's own data, and from their 'mastcrsct' differ by 
as much as 5 — 7 times the error propagated from the errors 
in photometry. It seems that there are systematic errors in 
at least one of these two data sets that have not yet been 
accounted for. Presumably this would affect any method to 
determine time delays and so one should not be surprised to 
find a spread in time delay results in the literature that is 
larger than the formal errors quoted. 
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APPENDIX 
SERIES 



CONTAMINATION OF THE TIME 



The analysis in the main paper is based on the assump- 
tion that the two time series can be measured without any 
contamination by foreground sources or instrumental off- 
sets or drifts, so that only random (measurement) noise is 
a source of errors. If one or both images of the quasar are 
contaminated by a foreground source, for instance the lens- 
ing object itself, this may influence the time lag determina- 
tion. If the foreground source is rapidly varying it can only 
give rise to false lags or aliases if that variability is signifi- 
cantly correlated with the time series of the lensed quasar. 
This is a priori unlikely because there is no causal physi- 
cal connection between the two objects. In practice it may 
occur if the measured time series is short compared to the 
time scales of variability of the lensed quasar and the hy- 
pothetical contaminating foreground source or if the phys- 
ical mechanisms causing the variability in the lens and the 
quasar have matching characteristic time scales. In the ab- 
sence of such correlation the variable part merely behaves 
as an extra somce of noise in the time series which can in- 
crease the uncertainty in the determined lag, but it will not 
have a systematic effect. A constant or low frequency non- 
zero contamination in one or both images can influence the 
determination of the time lag indirectly. It is demonstrated 
in this appendix that this can be detected if it occurs and 
can be corrected for. 

Consider again equations (2). Now, instead of having 
measured and F^''\ contaminated time series F^") and 
F^'' arc measured : 



F 



(a) 



(Al) 



m which C^"^ and C^*"' are the (unknown) constant contam- 
inating source (s) and the terms with coefficients L represent 
a possible (linear) drift. With these contaminated time series 
the kernels are built, and the inversion is carried out. In the 
main paper the normalization of the base functions is taken 
implicit for convenience of notation. Here these factors are 
be kept separate as coefficients {ui} or, for the case of the 
contaminated series, {ui} : 



'(b) 



Tmax \ 

J F(")(ti+T)dT 



(A2) 
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These normalization factors are calculated at the start of the 
algorithm and hence are known. For convenience of notation 
the following linear combinations g and h are defined : 
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-h) 





(A3) 



(A4) 



For the zero-order moment the following holds : 



(A5) 



Combining (A5) with (Al) yields : 

E-^Oa)~ia) pia) ^ ^ _ Ma) (Oa) _ J^ia) lAOa) 
^ I max 

^~*0a)~(a)^(6) ^ ^-^Oa)~(a) p(b) _^ (j{b) ^{Oa) (A6) 

Working through the equivalent of equations (8) it can be 
seen that instead of determining the constant / what is de- 
termined with the contaminated series is : 

1 



7(6) " I +^^--^5 



\ ^ I max"' 



(A7) 



1 /■(i') 
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(A8) 



For the linear target kernel T = t the coefficients c^^"^ 
satisfy : 

Y^-^W = r 

T^^' = 
Combining (A8) with (Al) yields : 

The same holds of course when the role of the time series is 
reversed. The resulting first order moments are : 

/to'") = Ito + ff^'"' [C*") - JC'")] 



to to 

7m 



j-g 



(16) 



fj(a) _ £^(6) 



(AlO) 



(la) 



Clearly a systematic deviation is present due to the contam- 
inating source in both the zero order and first order mo- 
ments. The presence of noise in the measured time series 
means that the equality in the first of equations (A5) and 
of (A8) become approximate just as in equations (7). This 
does not affect the determination in any systematic way. 

The equations (A7) and (AlO) form a system of four 
equations in the four unknowns I, to, yi = C'") — C^^^ /I, 
and 2/2 = L^'^'' — L^*''^ /I. Note that taking y\ and j/2 in this 
way refiects that with this method it is impossible to tell 
which of the two time series are contaminated, or to what 
extent either of them is contaminated. 

Since this system of equations is non-linear it is possi- 
ble that no solution exists or more than one solution. In any 
case the propagation of the data errors could cause large un- 
certainties in the parameters. Therefore the procedure that 
is followed is to solve (A7) and (AlO) to determine yi and 
2/2. The appropriate contamination is then subtracted from 
one of the two time scries and the inversion procedure is 
carried out again. For data free of random errors this sec- 
ond iteration should be contamination-free. For real data a 
few iterations can be necessary to reduce the contamination 
to zero within the measurement errors and an accurate de- 
termination of the relative magnification factor I and the 
time delay to is then possible directly from the inversion as 
outlined in the main part of this paper. 
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